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This  grant  provided  support  for  a  postdoc,  Lyudmyla  Barannyk,  at  the  University  of 
Michigan,  during  the  period  September  1,  2006  -  May  31,  2007,  to  assist  in  developing  a 
grid-free  particle  method  for  electrostatic  plasma  simulations. 


1  Objectives 

This  project  is  part  of  a  larger  effort  funded  by  AFOSR  (FA 9550-05- 1-0199,  Major  David 
Byers),  in  collaboration  with  Andrew  Christlieb  at  Michigan  State  University,  that  aims 
to  develop  a  grid-free  particle  method  for  plasma  simulations.  The  majority  of  plasma 
simulations  currently  use  the  particle-in-cell  (PIC)  method  [1,2],  The  advantage  of  PIC  is 
that  it  uses  a  fast  Poisson  solver  to  obtain  the  electric  field,  but  the  results  may  exhibit 
mesh-induced  effects  such  as  artificial  diffusion  and  anisotropy.  In  addition,  geometrically 
complex  domains  present  a  problem  for  the  fast  Poisson  solvers  typically  used  in  PIC* 

The  present  work  seeks  to  overcome  these  difficulties  and  substantially  improve  the  ac¬ 
curacy  and  efficiency  of  plasma  simulations  for  a  wide  range  of  applications.  Our  approach 
is  based  on  the  Lagrangian  formulation  of  charged  particle  dynamics,  in  contrast  to  the 
standard  Eulerian  formulation  in  terms  of  the  Vlasov- Poisson  equation.  The  Lagrangian  for¬ 
mulation  leads  naturally  to  a  grid- free  particle  method.  We  incorporate  several  techniques 
from  the  study  of  vortex  sheet  motion  in  computational  fluid  dynamics  [3]. 


2  Status  of  Effort 

The  investigators  made  progress  in  developing  the  Lagrangian  formulation  of  charged  par¬ 
ticle  dynamics,  deriving  a  regularized  system  that  satisfies  charge  neutrality  and  periodic 
boundary  conditions,  implementing  an  adaptive  particle  insertion  scheme,  and  performing 
test  computations.  The  work  is  being  written  up  for  publication  [4],  The  following  section 
summarizes  the  approach  and  presents  numerical  results  showing  the  method’s  capability. 
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3  Accomplishments/New  Findings 

Here  we  summarize  the  accomplishments  and  new  findings  obtained  during  the  award  period. 
Complete  details  are  in  an  article  being  prepared  for  publication  [4] .  We  start  by  recalling  the 
standard  Eulerian  formulation  of  charged  particle  dynamics  in  terms  of  the  Vlasov- Poisson 
equations. 

3.1  Eulerian  Formulation 

Let  f(x,  v,  t)  be  the  electron  probability  density  function  (pdf)  in  phase  space,  satisfying 

/(*.  v,  t)  >  0,  f  [  f(x,v,t)dxdv  =  1,  (1) 

J  — oo  j0 

where  x  is  space  (here  ID),  v  is  velocity,  and  t  is  time.  The  pdf  evolves  by  the  Vlasov 
equation, 

ft  +  vfx  -  Efv  =  o,  (2) 

where  the  electric  field  E(x,  t)  is  given  by  E  =  -<px  and  the  potential  function  d>(x,  t)  satisfies 
the  Poisson  equation, 

-<Pxx  =  P,  (3) 

with  periodic  boundary  conditions.  The  charge  density  is 

p(x,t)  =  -  [  f(x,v,t)dv  +  l.  (4) 

The  Vlasov-Poisson  equations  (2-4)  describe  the  evolution  of  the  pdf  f(x,  v,  t).  Next  we 
describe  an  equivalent  Lagrangian  formulation. 


3,2  Lagrangian  Formulation 

Consider  the  charge  flow  map, 

(a\  fx{a,P,t)\ 

\0J  \v(a,/M)/  ’ 


(5) 


where  (o,  0)  are  Lagrangian  parameters  labeling  a  charged  particle  at  location  (x,  v)  in  phase 
space  at  time  t .  This  is  shown  schematically  in  Figure  1 . 

The  equations  defining  the  flow  map  are, 


xt(a,0,t)  =  v{a,0,t),  (6) 

vt(a,0,t)  =  —f  f  (k(x(at0,t)^x{a,0,t))  +  x{a>p,t))u>a(a,0)dadl3  +  x(Q,0,t),  (7) 

J -oo  JO 


where 


*(x,  y)  =  —  sign(m  -  y ) 


(8) 


is  t.hc  electric  field  kernel  and  wo(a,  P)  is  the  initial  charge  distribution,  Equations  (6-7)  are 
simply  Newton’s  equations  in  first  order  form.  The  integral  on  the  right  side  oi  (7)  gives  the 
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Figure  1:  Charge  flow  map,  (a,/?):  parameter  space,  (or, u):  phase  space. 

electric  field  in  terms  of  the  flow  map  and  initial  charge  distribution.  We  assume  here  that 
the  spatial  domain  is  ID  and  the  solution  is  periodic  with  period  L  =  1,  but  these  are  not 
essential  restrictions.  Equations  (6-8)  are  the  desired  Lagrangian  formulation  for  charged 
particle  dynamics  in  phase  space. 

3.3  Particle  Discretization 

We  consider  a  cold  stream  in  which  the  charge  is  supported  on  a  curve  in  phase  space,  so 
that  the  flow  map  reduces  to  the  form  x(a,  t),v(a,  t).  We  use  the  midpoint,  rule  to  discretize 
the  electric  field  integral  in  (7).  Then  at  =  (i  -  \)/N,i  =  1  :  JV  are  the  discretization  points 
in  parameter  space  and  (.Tj{t),Ui(())  =  (ar(ai,t),u(ai,t)}  are  the  corresponding  particles  in 
phase  space.  The  discrete  evolution  equations  are 

Vi,  (9) 

N 

-  +  Xj)v)j  +  Xi,  (10) 

3= 1 

where  wj  =  l/N  are  the  quadrature  weights.  The  discontinuity  in  the  kernel  is  handled  by 
setting  k(x,  x)  =  0.  The  4th  order  Runge-Kutta  method  is  used  for  timestepping. 

Figure  2  shows  the  numerical  results  for  a  perturbed  uniform  stream  with  N  =  50,  At  = 
0.004.  The  solution  is  plotted  over  two  periods,  0  <  x  <  2.  Two  types  of  particles  are  plotted, 
active  (red)  and  passive  (blue).  The  active  particles,  (x;,i>j),  carry  charge  and  contribute  to 
the  sum  defining  the  electric  field.  The  passive  particles  have  Lagrangian  parameter  values 
at  the  endpoints  of  the  intervals  in  parameter  space.  The  passive  particles  are  used  in  the 
adaptive  particle  insertion  scheme,  described  below. 

The  stream  rolls  up  into  a  vortex  in  phase  space.  The  vortex  is  associated  with  particle 
trapping  as  particles  cycle  back  and  forth  in  the  given  period.  Other  particles  escape  and 
are  adverted  into  neighboring  periods,  leading  to  the  formation  of  long  thin  filaments  in  the 
stream.  The  simulation  loses  accuracy  as  time  proceeds,  as  shown  by  the  self- intersection 
of  the  stream  at  late  times.  Simply  using  a  smaller  time  step  A(  and  a  larger  number  of 
particles  /V  is  not  a  practical  way  of  maintaining  resolution.  Moreover,  as  shown  in  Figure  3a, 
even  as  the  computation  is  refined,  the  inner  portion  of  the  spiral  is  tangled.  We  believe 
this  tangling  is  due  to  a  singularity  in  the  exact  solution.  To  address  these  issues  we  employ 
regularization  and  adaptive  particle  insertion,  as  explained  in  the  next  two  sections. 
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Figure  2:  Cold  stream,  N  =  50,  At  —  0.004,  active  particles  (red),  passive  particles  (blue). 
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Figure  3:  Cold  stream,  closeup  of  spiral  core  at  t  =  8,  (a)  The  solution  converges  slowly 
as  N  increases  and  At  decreases.  The  inner  portion  of  the  spiral  remains  tangled,  (b)  The 
regularized  solutions  with  5  >  0  roll  up  smoothly  and  are  free  of  tangling. 
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3.4  Regularization 

Consider  the  following  regularized  approximation  of  the  electric  field  kernel, 

k‘{z’v)  “  ?((x-»r~+ *»)■/»'  (11) 

where  S  >  0  is  a  smoothing  parameter.  The  exact  kernel,  in  equation  (8),  is  recovered  in  the 
limit  5  — *  0.  Using  the  regularized  kernel  we  obtain  a  regularized  system, 

x({a,  t)  =  u(a,t),  (12) 

/DO  fl  r  1 

/  fca(x(Q,C),x{Q,t))wo(d)dd  +  p  /  ks(x(a ,t),y)dy  -  a,  (13) 
DO  Jo  J  0 

where  the  constants  p  and  a  are  chosen  to  enforce  charge  neutrality  and  periodicity  [4],  Using 
the  midpoint  rule  as  before,  we  obtain  a  particle  discretization  of  the  regularized  system. 
As  shown  in  Figure  3b,  the  solution  of  the  regularized  system  is  free  of  tangling  for  S  >  0. 
But  there  is  still  a  need  for  adaptive  particle  insertion  to  maintain  resolution  without  using 
prohibitively  large  N. 

3.5  Adaptive  Particle  Insertion 

To  maintain  resolution  as  the  curve  evolves,  we  adaptively  insert  new  particles  at  each  time 
step  [3],  Figure  4  defines  two  quantities,  (a)  di.  chord  length  of  the  interval,  (b)  d2:  distance 
from  the  active  particle  to  the  chord.  If  either  of  these  quantities  is  larger  than  a  user-specific 
value,  the  interval  is  split,  as  depicted  in  Figure  4c.  After  an  interval  is  split,  the  quadrature 
weights  are  reset. 


Figure  4:  Adaptive  particle  insertion,  (a)  d^:  chord  length,  (b)  d2:  particle-chord  distance, 
(c)  panel  splitting;  the  active  particle  becomes  passive  and  two  new  passive  particles  are 
inserted  using  quadratic  interpolation  with  respect  to  the  Lagrangian  parameter. 


6 


3.6  Numerical  Results 

Figure  5a  shows  a  simulation  of  the  regularized  system  (11-13)  with  <5  =  0.1  and  adaptive 
particle  insertion.  The  initial  number  of  particles  is  N  =  400  and  the  final  number  is 
N  =  2920.  To  clarify  the  charge  transport  in  phase  space,  the  charge  coming  from  each 
initial  period  is  colored  (e.g,  blue:  0  <  x  <  1,  red:  1  <  x  <  2,  etc.).  The  closeup  in 
Figure  5b  shows  that  the  fine  scale  structures  are  well-resolved  and  free  of  tangling.  The 
results  suggest  the  presence  of  chaotic  dynamics  due  to  a  heteroclinic  tangle. 

Figure  6  shows  a  preliminary  simulation  of  the  cold  two-stream  instability.  In  this  case 
there  are  two  streams  of  particles  moving  in  opposite  directions.  The  initial  charge  density 
was  perturbed.  The  streams  roll  up  smoothly  into  a  vortex  in  phase  space.  This  is  a  step 
towards  extending  the  code  to  handle  warm  charge  distributions. 

3.7  Ongoing  Work 

(a)  We’re  performing  various  checks  of  the  numerical  accuracy,  e.g.  energy  conservation,  time 
integration  (for  5  =  0  and  6  >  0). 

(b)  We  believe  that  a  cusp  singularity  forms  at  a  finite  time  in  the  charge  distribution  and 
electric  field  as  a  function  of  the  spatial  ^-coordinate.  We  will  document  this  numerically 
using  spectral  analysis  (FFT). 

(c)  We’re  studying  the  phenomenon  of  particle  trapping.  The  aim  is  to  determine  how  much 
of  the  charge  remains  trapped  in  its  original  period,  how  much  is  advected  and  trapped 
in  neighboring  periods,  and  how  much  undergoes  free  streaming.  Charge  transport  is  a 
fundamental  issue  and  the  present  approach  rnay  well  provide  new  insights. 

(d)  The  code  is  being  extended  to  handle  warm  distributions,  i.e.  those  in  which  the  charge 
is  distributed  over  a  region  in  phase  space,  as  opposed  to  simply  a  curve  or  a  finite  set  of 
curves.  This  will  enable  us  to  study  Landau  damping  and  compare  with  previous  results. 

(e)  We  intend  to  perform  a  detailed  comparison  with  PIC  simulations. 


3.8  Summary 

The  results  presented  in  the  previous  sections  demonstrate  the  capability  of  the  grid- free 
particle  approach  under  development.  The  method  is  especially  well  suited  for  studying 
charge  transport  and  lamentation.  We  believe  this  approach  will  substantially  improve  the 
accuracy  and  efficiency  of  plasma  simulations  for  a  wide  range  of  applications. 


4  Personnel  Supported 

The  project  provided  50%  support  for  a  postdoc,  Lyudmyla  Barannyk.  Dr.  Barannyk  is 
starting  a  tenure-track  position  at  the  University  of  Idaho  in  August,  2007.  The  project 
extends  the  work  started  last  summer  by  Benjamin  Sonday  when  he  was  an  undergraduate. 
Mr.  Sonday  is  currently  a  graduate  student  at  Princeton  University  with  support  from  a 
Department  of  Energy  Computational  Science  Graduate  Fellowship. 
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Figure  5:  Cold  stream,  long  time  evolution,  regularized  (S  =  0.1),  adaptive  particle  insertion. 
Charge  coming  from  each  initial  period  is  colored,  e.g.  blue:  0  <  x  <  1,  red:  1  <  x  <  2,  etc. 
(a)  time:  0  <  t  <  20,  initial  N  =  400,  final  N  =  2920,  (b)  closeup  at  t  =  20, 
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Figure  6;  Cold  two-stream  instability. 


5  Publications 

The  work  accomplished  during  the  support  period  is  being  prepared  for  publication  [4],  The 
article  will  be  sent  to  the  program  manager,  Dr,  Fariba  Fahroo,  upon  completion.  The  group 
has  published  several  previous  articles  with  AF  support  [5-7]. 

6  Interactions/Transitions 

The  work  supported  by  this  award  was  presented  at  three  conferences. 

•  American  Physical  Society,  Division  of  Plasma  Physics,  Philadelphia,  November,  2006 

•  Society  for  Industrial  and  Applied  Mathematics,  Conference  on  Computational  Science 
and  Engineering,  Costa  Mesa,  February,  2007 

•  International  Congress  on  Industrial  and  Applied  Mathematics,  Zurich,  July,  2007 

It  will  also  be  presented  at  the  upcoming  APS  meeting  of  the  Division  of  Plasma  Physics, 
Orlando,  November,  2007.  The  project  is  part  of  a  larger  effort  supported  by  AFOSR 
(FA9550-05-1-0199,  Major  David  Byers)  in  collaboration  with  Andrew  Christlieb  at  Michigan 
State  University. 
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